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1.0 INTRODUCTION 


1.1 General 

VSAERO is a computer program for calculating the subsonic 
AEROdynamic characteristics of arbitrary configurations having 
Vortex Separation and strong Vortex/Surface interaction. This 
document describes the theory and numerical procedures used for 
the method. An earlier report (1) describes the input and use of 
the computer program. 


1.2 Objectives 

The initial objectives of the new program were to analyze 
general wing configurations with multiple part-span, high-lift 
devices and to give special attention to edge-vortex separations 
and to close vortex/surface interactions. Later objectives 
brought in the modeling of fuselage and canard surfaces leading 
to the present capability for arbitrary configurations including 
complete aircraft. 


1.3 Strategy 

At the beginning of the program development there was only 
one clear cut choice for approaching the solution to the above 
nonlinear problem. This was an iterative viscous/potential flow 
calculation using a potential flow panel method coupled with 
integral boundary layer methods, with wake-shape iteration calcu- 
lations included in the potential flow analysis. No other ap- 
proach appeared to offer the feasibility of achieving the above 
objectives. Finite difference field methods, for example, 
treating more general flow equations, could not be developed 
within a reasonable time scale for application to the complex 
geometry of multiple part-span, high-lift devices. In the mean- 
time, the surface integral approach of the panel method offered a 
powerful and practical engineering tool which had the capability 
to represent the complex geometries involved in the objectives. 
Furthermore, the low computing time of the panel method was an 
attractive feature in view of the complex iterative loops made 
necessary by the nonlinear effects of viscosity and wake loca- 
tion. 


In the chosen strategy for the development of the VSAERO 
program the capabilities of two earlier programs were combined: 
(i) VIP3D (2) offered a viscous/potential flow iterative method 
for transport wings with full-span, high-lift devices; and (ii) 
QVORT, a quadrilateral-vortex panel method (3) and (4), offered 
an iterative wake relaxation procedure. Neither method had a 
completely satisfactory singularity model or geometry package for 
the new objectives. Accordingly, in the development of VSAERO 
these features were given prime consideration. 
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Various singularity models were investigated for the new 

F rogram including symmetrical singularities (5), high-order 
ormulations and numerical integration. Finally, a low-order 
formulation based on internal Dirichlet boundary conditions was 
chosen. This formulation is described in Section 2. The numeri- 
cal procedure for the method, which employs quadrilateral panels 
of constant doublet and source distributions, is described in 
Section 3. The description includes details of the extended 
geometry package for the VSAERO program. Correlation studies are 
examined in Section 4; these include a number of basic test cases 
involving exact, nominally exact and also experimental data. 
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2.0 FORMULATION 

2.1 Flow Equation 

In order to obtain a solution to the flow problem the sur- 
face of the configuration is assumed to be at rest in a moving 
sirstream. Regions of the flow field that are dominated by 
viscous and rotational effects are assumed to be confined to thin 
surface boundary layers and wakes. The rest of the flow field is 
assumed to be inviscid, irrotational and also incompressible. 

Figure 1 shows a streamwise section taken through a wing and 
its wake. The section indicates the surface, S, of the wing, the 
surface (upper and lower) of the wake, W, and an outer surface. 
So,, which encloses the complete flow problem at infinity. The 
total surface, S + W, divides space into two regions: the exter- 
nal region which contains the flow field of interest and the 
inner region which contains a fictitious flow. Velocity poten- 
tial fields, $ and , are assumed to exist satisfying Laplace's 
equation in the flow field and in the inner region, respectively; 
i.e. , 


and 


V 2 <J> = 0 in the flow field 


V 2 $ . 

l 


= 0 


in the inner region. 


Next, apply Green's Theorem to the inner and outer regions 
and combine the resulting expression (see Lamb (6)). This yields 
the following expression for the velocity potential, $p, anywhere 
in the two regions, expressed in terms of surface integrals of 
the velocity potential and its normal derivative over the bound- 
ary surface: 


St- WH S 


* ± ) n • V(^)dS 



st-wt-s 


(V$ - v$ i ) ds 


( 1 ) 


where r is the distance from the point, P, to the element, dS, on 
the surface and n is the unit normal vector to the surface 
pointing into the fluid. 


The first integral in Eq. (1) represents the disturbance 
potential from a surface distribution of doublets with density ($ 
- $j_) per unit area and the second integral represents the con- 
tribution from a surface distribution of sources with density 
-n»(v$ - per unit area. The singularities, therefore, rep- 
resent the jump in conditions across the boundary; the doublet 
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density represents the local jump in potential and the source 
density represents the local jump in the normal component of 
velocity. 


Examine Eq. (1) separately over the three surfaces, S, W, 
S m , as indicated in Figure 1. 

The surface at infinity, S^, can be regarded as a large 
sphere centered on P. The local conditions at that boundary 
consist only of the uniform onset velocity, V ro , the disturbance 
at that distance due to the configuration having essentially 
disappeared. The surface integrals in Eq. (1) taken over then 
reduce to 4>oo p ; i.e., the velocity potential at point, P, due to 
the onset flow. 

The upper and lower wake surfaces, W, are assumed to be 
infinitesimally close to each other for the present problem; 
i.e., thin wakes. (A simplified thick wake representation is 
described in References 7 and 8.) In this case the corresponding 
upper and lower elements can be combined and the 4>j_ term for this 
part of the integral disappears. Furthermore, if we neglect 
entrainment into the wake surface, then the jump in normal 
velocity component, n • (V<J>u - V$ L ) , is zero and so the source 
term in Eq. (1) disappears for the wake contribution. (Note that 
entrainment modeling could be included given a distribution of 
the entrainment normal velocity.) 


With the above conditions, Eq. 


(1) becomes: 


*P 


_ 1 _ 

4n 



$ i ) n ‘ 7 ( p )ds 



n • (V<|> - V<K )dS 


+ 



V (^)dW + d> 

L T oo 

P 


( 2 ) 


The normal, n, on surface, W, points upwards in this case. 

If the point, P, lies on the surface, S, then the integral 
becomes singular on S. To evaluate the local contribution the 
point, P, is excluded from the integral by a local spherical 
deformation of the surface centered on P. For a smooth surface 
the local deformation is a hemisphere and the local contribution 
obtained in the limit as the sphere radius goes to zero is 
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ij ($ - <j>^ ) p , which is half the local jump in potential across the 
surface at P. Note that the first potential in the expression is 
on the side of the surface on which P lies, or, more con- 
veniently, the sign can be changed for the point on the other 
side. For example, if P is on the inside surface of S, then Eq. 
(2) becomes: 


$ P 



S-P 


$ ± ) n • V(^)dS - h ( 3 > - $ . ) 

1 P 


+ 



(7$ - 7$ i )ds 


$ L ) n * V(^)dw + 


w 



(3) 


2.2 Boundary Conditions 

A solution to Eq. (3) must satisfy a number of known bound- 
ary conditions which can be imposed — or in the case of the wake — 
implied at the outset. On the surface, S, the external Neumann 
boundary condition must be satisfied. 

n • 7$ = -v N - n • V g on S (4) 

where V N is the resultant normal component of velocity relative 
to the surface. V N is zero for the case of a solid boundary with 
no transpiration; non-zero values are used to model boundary 
layer displacement effect, entrainment and also inflow/outflow 
for engine inlet/exhaust modeling. V g is the local velocity of 
the surface; this may be composed of several parts due to body 
rotation (e.g., pitch oscillation (9), or helicopter blade rota- 
tion (10)), arbitrary motion (11), body growth (12), etc. For 
the present problem the surface of the configuration is assumed 
steady in a body-fixed coordinate system and so V s is zero. 

Hence , 


n • 7$ = -V on S (5) 
N 

The wake surface, W, cannot support a load and so the doub- 
let distribution, $ - $ , on W must satisfy a zero-force condi- 
tion. u L 
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Consider an element, dW, of the wake surface with local unit 
normal, n. The local vorticity vector associated with the doub- 
let distribution, $ rT - 3> , has the density 

U Li 

y = -n x y ($ - ) (6) 

The elementary force exerted on the element, dW, in the 
presence of the local mean velocity, V, is (from the Kutta- 
Joukowski Law): 

6 F = pV x y dW 

where p is the local density. And since the force must be zero, 
it follows that 


V x y = 0. 

Substitute for y from Eq. (6) : 

V x (n x y (<j> - * )} = 0 

U L 


and expand 


n V • V($ - *_) - V • nV($ TT - 3> ) =0. (7) 

U Li U Ll 

Equation (7) is satisfied when V • n = 0; i.e., the surface, 
W, is aligned with the local flow direction, and V • 7(3^ - Or) = 
0. In the latter case, the gradient of the doublet distribution 
in the direction of the local mean flow is zero; in other words, 
the wake doublet distribution is constant along mean streamlines 
in the wake surface. The doublet value on each streamline is, 
therefore, determined by the condition at the point where the 
streamline "leaves" the surface, S. Clearly, at the outset of 
the problem with both the doublet distribution and wake location 
unknown, an iterative approach is necessary to obtain a solution. 
When a converged solution is obtained the upstream edge of the 
wake, W — and hence the trailing edge of the surface, S — carries 
no load and so the Kutta tr ailing-edge condition will be satis- 
fied. At the outset, therefore, the Kutta condition is im plied 
simply by shedding the trailing-edge potential jump (3^ - 3> L ) as 
a constant down each "streamwise" line on an initially prescribed 
wake surface. An explicit Kutta condition in which the doublet 
gradients on the upper and lower surfaces at the trailing edge 
are equated to cancel has been used on occasion, but in the 
formulations described below, this does not appear to be neces- 
sary. 
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2.3 Choice of Singularity Model 


In principle, a given flow field can be constructed from an 
infinite number of combinations of doublet and source distribu- 
tions on the surface, S, each combination producing a different 
flow in the inside region. In practice, however, there may be 
only a small number of singularity combinations which are suit- 
able as a basis for a well behaved numerical model; i.e., one 
that is accurate, convenient to use and robust (insensitive to 
user abuse). 

A unique combination of singularities can be obtained in Eq. 
(3) in a number of ways. One way is to specify one of the 
singularity distributions and to solve for the other using 
boundary conditions only on one side of the boundary. There are 
several examples of this: in the source-only formulation (13) 
the doublet distribution is set to zero; in the doublet-only 
formulation, e.g., (3) and (4), the source distribution is set to 
zero; in program VIP3D (2), the doublet distribution (or vortici- 
ty) is specified and the method solves for the source distribu- 
tion; the reverse of the latter has also been used, i.e., speci- 
fying the source distribution — usually related to the thickness 
distribution — and solving for the doublet or vorticity distribu- 
tion. 


Another way of achieving a unique combination of singulari- 
ties is to apply certain constraining relationships on the singu- 
larity distributions. The "symmetrical singularity method" (5) 
is one example of this and was examined briefly at the start of 
the VSAERO program development. 

One characteristic that separates "good" singularity models 
from "bad" ones from the numerical point of view is that the flow 
field generated in the inner region by a "good" model is rather 
benign and related to the boundary. This is usually not the case 
for a "bad" model. In other words, when passing through the 
boundary, S, the jump from the internal flow to the external flow 
should be small — thus requiring a minimum of perturbation from 
the singularities (Eq. (3)). "Bad" singularity methods were 
observed to have very large internal cross flow between the upper 
and lower surfaces of wings and required high panel densities in 
order to achieve acceptable accuracy in the flow solutions. This 
resulted in a push towards high-order formulation which, for 
subsonic flows at least, proved quite unnecessary. 

One way of achieving "good" characteristics is to treat the 
internal flow directly in Eq. (3). This is, in fact, another way 
of obtaining a unique singularity distribution — in this case by 
specifying boundary conditions on both sides of the surface, S. 
Earlier methods (14), (15) specified the velocity on the inside 
surface of S. Alternatively, the inner velocity potential, , 
can be specified directly in Eq. (3). This is referred to as an 
internal Dirichlet boundary condition. Three possible internal 
flows were examined for the VSAERO program. Figure 2 . Two of 
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these are described below, the third one, in which the internal 
flow was aligned with the wing chord line (e.g., = -xV^ cos ct ) 

was examined in a two-dimensional pilot program only. Although 
this worked very well, the formulation was not as convenient as 
the other two. Even so, it would be a useful member in a pos- 
sible general scheme in which the user could select a suitable 
internal flow for each component in a complicated configuration. 

(i) Zero Flow Inside 

Consider first the interior stagnation condition which was 
implied by earlier doublet-only codes (e.g.. References 3 and 4) 
and which used the exterior Neumann boundary condition, V • n = 
0. For the present formulation, set $ . = constant = 0, say in 

Eq. ( 3) , giving 1 

$ n • V(p)dS - 

S-P 




• V $ dS 




+ 



V 


n 


V(^)dW+ <j> c 


(8) 


The doublet distribution is now the total velocity potential 
at the surface and the source distribution is now the total 
normal velocity of the fluid at the surface, S. This must 

satisfy the external Neumann boundary condition (Eq. 5); i.e., 

n • V; = -v N 

This term can be used directly in the second integral of Eq. (8) . 

For the solid boundary problem, i.e., zero transpiration, V N 
is zero and the source term disappears leaving 



S-P 


n 


7(^)ds 


h<i> 


P 


+ 




- $ T > n 

-L 
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7(|)dW + <f> rop = 0 


(9) 
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This is by far the simplest formulation for the lifting case 
with no transpiration. It is, in fact, even simpler than the 
original zero-lift source method based on the external Neumann 
boundary condition (13). In a low-order formulation (i.e., con- 
stant doublets on flat panels) the method has some minor problems 
associated with numerical differencing for velocities and with 
the effect of panel arrangement on accuracy (16). These problems 
can be alleviated by separating $ into a known part and an un- 
known part — the latter being as small as possible. 

The most convenient breakdown is to use the onset flow 
potential as the known part; i.e., $ = (p^ + <f>. The known part, 

< Poo f and its velocity field, , are used directly in the formula- 
tion. The solution for <p and the subsequent numerical differen- 
tiation for surface velocity are then less prone to numerical 
error . 

For the general case with transpiration, inlet flows, jets 
and unsteady conditions, the source term must be present and so 
the advantage of the doublet-only formulation diminishes. In the 
meantime an alternative "zero internal perturbation" formulation 
had become popular in other methods and is considered below. 


(ii) Internal_Flow Equal to the Onset Flow 

The internal Dirichlet boundary condition, was used 
earlier by Johnson and Rubbert (17) and also Bristow and Grose 
(18) in high-order formulations. With <J>. = d> and = <b , Eq. 
(3) becomes: 1 p °° H 


<J> n • 7 (•£■) dS - h<p p + 

S-P 


£ n • (V* - V4> 00 )dS 

s 


where <f>, the perturbation potential on the exterior surface, is 
now the doublet density; i.e., 

4 TTy = tp = $ - ^ (11) 

The source distribution in this formulation is 

4 tto = - n • (V$ - V<j> ) 

' oo 




// 


($u - $ L ) n • V(-)dW 


w 


= 0 


( 10 ) 
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Expanding this expression and substituting for n • from the 
external Neumann boundary condition in Eq. (5): 

4no = V - n • v 

N 0 

The source distribution is therefore established at the outset. 


Again, for solid boundary conditions, Vjj is zero; for the 
more general case here, V N may have two parts representing (i) 
boundary layer displacement effect using the transpiration 
technique (2), and (ii) inflow/outflow for engine inlet/exhaust 
modeling . 


Thus , 


4 it a = (V 6*) 
9 s e 


+ V 


NORM 


- n 


( 12 ) 


The boundary layer term is set to zero for the first poten- 
tial flow solution; on subsequent iterations it is updated by 
coupled boundary layer calculations. The Vj, orm term has a posi- 
tive value for outflow and negative for inflow. 


Compared with the zero internal flow formulation (Eq. 9), 
the present formulation (Eq. (10) is more forgiving for "bad" 
panel arrangements (e.g.. Reference 16) and is the current formu- 
lation in the VSAERO program. The reason for its better behavior 
probably arises from the smaller jump in the flow condition from 
the inner to the outer flow; i.e., the singularity strengths are 
smaller. However, there are situations where this is not the 
case; e.g., a wing at large angle of attack or a powered nacelle 
in zero onset flow. In other words, it may be desirable to 
provide a capability that allows the internal flows to be speci- 
fied independently of the onset flow so that they can comply more 
closely with the shape of each component. 

Once the doublet solution is known, then Eq. (3) can be 
applied to determine the velocity potential field, the gradient 
of which gives the velocity field. For the present formulation, 
Eq. (3) becomes: 





• V(^)dS + Ky p 



n • V(-)dW + 4^ 

r p 


dS 


( 13 ) 
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Where the doublet, y, and source, o , are taken from Eqs. 
(11) and (12), respectively and uw = t * ie wa *^e doublet 
distribution. The factor, K, has different values depending on 
the location of the point, P; if P lies off the surface, then K = 
0; if P lies on a smooth part of the surface, then K = 2tt if p is 
on the outside, and K = -2tt if P is on the inside; if P is at a 
crease in the surface, then K takes the value of the appropriate 
solid angle contained at the crease. In each case when P is on 
the surface, then the surface integral terms in Eq. (13) exclude 
P. 


If the potential field has been computed at a mesh of points 
then the velocity field can be generated using local numerical 
differentiation; i.e.. 


Vp = - v$p (14) 

Alternatively, the velocity field can be obtained by taking 
the gradient of Eq. (13) directly with respect to the coordinates 
of point P, thus forming the source and doublet velocity kernels. 
(This is the present approach in VSAERO.) 



w 


o 7 Op) dS 
S 

)dW + V^, (15) 



(The kernels will be developed further in Section 3.) The off- 
body velocity field is used in the wake-relaxation cycle, in off- 
body streamline calculations and for general flow-field informa- 
tion. 


The discussion so far has been concerned with thick con- 
figurations having a distinct internal volume enclosed by the 
surface, S. If parts of the configuration are extremely thin 
(e.g., thickness/chord ratio < 1%, say) or are wing-like and 
remote from the area of interest, then these parts may be repre- 
sented by open surfaces. When the upper and lower parts of the 
surface are brought together, the corresponding upper and lower 
elements can be combined and Eq. (2) becomes 
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*p 


4 TT 




V 


n 


V (^)dS 



n • (V$ 0 - 7$ l ) dS 


+ 


4 7T 



($U - * L ) 


n 


V (^) dW + cf> ro 


(16) 


If the normal velocity is continuous through the sheet then 
the term n • (V'I’u ~ V$l ) = 0 and the source term disappears. 
(This is act a necessary step; the source term could be left in 
the equation for simple modeling of thickness effects.) Removing 
the source term does not preclude the use of non-zero normal 
velocity at the surface; the only restriction is that the normal 
velocity be continuous through the surface. Equation (16) can 
then be written: 



V (^)dS 




7 (^)dW + ^ (17) 

P 


where y = ^ - # L ) is the jump in total potential across the 

sheet . 


In order to satisfy the external Neumann boundary condition 
(Eq. (5) ) the velocity expression is required. Hence, taking the 
gradient of Eq. (17) with respect to the coordinates of P, 


V P= " 



y V (n • V (~) )dS 


I ( % v<n • * ^) )dW + V, 


(18) 


W 


And, applying the boundary condition, = V„ from Eq. (5) 

then p p N p 
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• V (n 


• V(I))dS 
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W 


n 


V(n 


V ( -p) ) dW + n p 



(19) 


This is the basic equation for the unknown doublet distribu- 
tion on thin surfaces (referred to as Neumann surfaces in VSAERO 
in accordance with their primary boundary condition). Again, the 
velocity kernels will be developed further in Section 3. 

In the VSAERO program, thin and thick components may exist 
simultaneously in a configuration. 
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3.0 


NUMERICAL PROCEDURE 


The treatment of a given flow problem can be broken down 
into a number of distinct steps, each step having a certain 
function. These steps are shown in the VSAERO program outlined 
in Figure 3. The first step is to define the surface geometry of 
the configuration and to form the panel model. A m atrix of 
influence coefficients is then formed; this represents a set of 
simultaneous linear algebraic equations resulting from a discre- 
tization of the flow equations in Section 2.0. The solution to 
these equations provides the surface doublet distribution from 
which the surface velocities and, hence, surface streamlines, 
pressures and body forces and moments can be analyzed . With the 
surface singularity distributions known, the wake-shape calcula- 
tion can proceed, using computed off-body velocities. The new 
wake shape must then be repanelled and the matrix of influence 
coefficients modified for the next pass through the wake-shape 
iteration loop. 

The surface flow solution can be passed to the boundary 
layer calculation. The computed boundary layer displacement 
growth rate is passed back to the matrix procedure where the 
right-hand side terms of the equation are modified for the next 
pass in the viscous/potential flow iteration loop, Figure 3. 
Finally, when both iteration loops have been completed, off-body 
surveys for velocity and streamline paths are processed. 

The flow chart showing the relationships of subroutines 
which perform the above steps is given in Appendix A. Details of 
the steps are described in the following Subsections. 


3.1 Geometrical Description 

The surface geometry of the configuration to be analyzed is 
described in a global coordinate system (G.C.S.) which is a body- 
fixed frame of reference, Figure 4. The configuration is divided 
into a number of parts for convenience. Figure 5. The functions 
and use of these parts are described in detail in the Program 
User's Manual (1). Briefly, the two major parts are the (solid) 
SURFACE and the (flexible) WAKE. On a complicated configuration 
there may be several surfaces and several wakes. The surface and 
wake are further subdivided as described below. 


3.1.1 Patch 

Each surface is represented by a set of quadrilateral PANELS 
which are assembled into a regular array of rows and columns 
within PATCHES, Figure 6. Patches are grouped together under two 
headings (Figure 5), COMPONENTS and ASSEMBLIES for user con- 
venience (1). 
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Figure 3. Outline of Steps in the VSAERO Program. 
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Figure 5. Configuration Hierarchy. 
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PANEL SEQUENCE PANEL COLUMN 
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Figure 6. Patch Conventions. 


Each patch has a user-supplied IDENT parameter. IDENT=1 and 
2 identify patches on wings and fuselages, respectively; the main 
difference at this time is the printout of computed 
characteristics (see Subsection 3.4). Both of the above IDENT 
values are associated with thick surfaces and use the internal 
Dirichlet boundary condition. Patches with IDENT=3 are asso- 
ciated with thin (i.e., open) surfaces on which only the Neumann 
boundary condition is used. Again, the main difference as far as 
the user is concerned is a different printout of computed charac- 
teristics, in this case, pressure and velocity information is 
provided on both sides of these surfaces. 

A patch is basically a four-sided shape when its surface is 
developed; i.e., "opened out". Because the initial objectives 
were directed towards wing surfaces, the terms "chordwise" and 
"spanwise" were adopted to describe the two directions across the 
patch. Figure 6. These terms are still used within the patch 
even when the patch has a more arbitrary orientation, e.g., on a 
fuselage. 

The order of the sides of a patch must proceed in the anti- 
clockwise direction when viewed from the flow field. Figure 6. 
This is important as it affects the direction of the computed 
surface normal vector which must point into the external flow 
field (Subsection 3.1.2). 

Patch geometry is defined using a set of "chordwise" lines 
called SECTIONS. The set of sections distributed in the "span- 
wise" direction across the patch define its surface. Figure 7. 
The points that are input to describe a section geometry are 
referred to as BASIC POINTS and are not necessarily panel corner 
points. From these points the program generates by interpola- 
tions a set of TEMPORARY CHORDWISE POINTS using information 
provided by the user on NODE CARDS (1). The actual panel corner 
points are obtained by interpolation along spanwise lines; these 
lines pass through a corresponding temporary chordwise point on 
each section. The interpolation scheme used in both directions 
uses a biquadratic equation and is described in Appendix B. The 
interpolation scheme is part of an automatic paneling procedure, 
use of which is optional: the user has the choice to input the 
actual panel corner points if he wants to. A number of input 
options are provided in the program to help the user generate the 
panel model (1). 


3.1.2 Panel Geometry 

The four corner points, R^, i=l,4, defining a quadrilateral 
panel are in the same sequence as the corners of the parent 

patch. Figure 8. Two diagonal vectors are constructed from these 
points : 
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SPANWISE 

SEQUENCE OF BASIC POINTS 



Figure 7 . 


Sections Defining Patch Surface 



Di — R3 — 


D 2 R 4 R 2 

The vector product of these diagonals produces a vector 
normal to the mean plane of the quadrilateral; i.e., 

n = D 2 x d 2 / ID x x D 2 I 

The order of the sides, therefore, is important to ensure 
that the unit normal is always directed outwards from the surface 
into the flow field. The modulus of the diagonal vector product 
also provides the area of the quadrilateral projected onto the 
mean plane: 

AREA = IDj x d 2 I /2 

The center point is defined as the mean of the four corner 
points : 

r c= 1 J2 R i 

i=l 

This point is also used as the panel's control point where 
the boundary conditions are satisfied (see Subsection 3.2). 

Two unit tangent vectors, £, m, are constructed, which, 
together with n form a right-handed orthogonal unit vector system 
for local coordinates, Figure 8. The origin of this system is 
the center point, R c . 

Tangent vector, m, is directed from R c to the mid-point of 
side 3 of the quadrilateral: 

m = { (R + R 4 >/2 - R c >/ I (R 3 + R 4 >/2 - R c l 

Thus id lies in the mean plane of the quadrilateral even if 
the corner points are not co-planar. 

The tangent vector, £ , is constructed orthogonal to m and n: 

£ = m x n 

The four corner points are projected onto the mean plane and 
expressed in terms of the local coordinate system. For example, 
the first point becomes 

E 1 = R 1 “ R c 
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which has components, E r , E , in the i and m directions, 
respectively, where 1 1 



The vertical offset, • n, is the projection distance of 
the corner points from the mean plane and is a measure of the 
skew of the quadrilateral. On a skewed quadrilateral the four 
corner points are equidistant from the mean plane, two being 
above and two below. The magnitude of this offset should be kept 
small in relation to the size of the panel to avoid large holes 
in the singularity panel representation of the surface. The 
projected (flat) panel is used in the evaluation of influence 
coefficients (Subsection 3.2). 

Finally, the half median lengths, SMP, SMQ, are evaluated 
for the quadrilateral. These are (Figure 8) 

SMP = I (R 0 + RJ/2 - R I 

Z 3 c 


SMQ = I (R^ + R 4 )/2 - R c l 


These are the half-lengths of the diagonals of the paral- 
lelogram which is always formed when the mid-points of the 
adjacent sides of a quadrilateral are joined (even if the lat- 
ter's corner points are not co-planar); the parallelogram lies in 
the mean plane of the quadrilateral and its area is half that of 
the projected quadrilateral. 

Within each patch, the regular arrangement of panels causes 
the adjacent side mid-points of neighboring panels to coincide 
exactly. This allows the SMP and SMQ lengths of panels to be 
linked, respectively, in the chordwise and spanwise direction 
over the patch and thereby provides a close approximation to 
surface distances between panel centers. This feature is useful 
in the evaluation of the perturbation velocities (i.e., the 
derivative of u with respect to surface distance, see Subsection 
3.4) . 

Finally, the quantity (SMP + SMQ) is used as the average 
"size" of the panel when expressing the distance of a point from 
the panel in terms of panel size. 


3.1.3 Panel Neighbors 

In order that a reasonable two-way differentiation of the 
surface doublet distribution can be performed (see Subsection 
3.4), it is important to locate quickly for each panel the set of 
four neighboring panels and their orientation. Each panel, 
therefore, keeps an array of four neighboring panels, NABOR^, 
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i-1,4 (in the same sequence as its sides), together with the 
adjacent side numbers of those panels, NABSID^, i=l,4. Figure 9. 

Clearly, within the regular grid system of a patch, locating 
neighbors is easy; even so, the neighbor information is still 
stored to form a consistent system and to avoid repetitive calcu- 
lations. Also, this allows neighbor relationships to be broken 
between pairs of panels, when, for example, a wake is shed from 
their common edge. This prevents a doublet gradient being 
evaluated across a line where there is an actual jump in poten- 
tial (see Subsection 3.4). 

Across the joints between patches panel neighbors are not 
immediately available; for example, one panel may be neighbor to 
several smaller panels on an adjoining patch. In order to locate 
these neighbors an automatic procedure has been installed in the 
code which scans patch side panels in a search for possible 
neighbors. In this search only patches within the same assembly 
are considered. "Undesirable" neighbors are quickly eliminated 
on the basis of relative geometry during the assembly of a short 
list of possible neighbors for each side panel. From this list 
of candidates one PREFERRED NEIGHBOR is selected. At this time, 
preference is given to the panel whose control point lies closest 
to a normal plane constructed on the side panel. Figure 10. This 
plane contains the side panel's control point, unit normal vector 
and the side mid-point at the patch edge. 


3.1.4 Wakes 

A wake consists of a number of quadrilateral wake panels 
which are assembled in streamwise columns. Each column is as- 
sociated with an upper and a lower wake shedding panel with 
subscripts, KWPU, KWPL, respectively. Figure 11. The common edge 
of these two panels is the upstream edge of the wake column, the 
VSAERO program coding allows any panel edge in the configuration 
to shed a wake (l)--it is up to the user to specify a wake model 
that is physically realistic. 

The streamwise edges of each column are identified as WAKE 
LINES, Figure 12. These lines are specified by the user to 
describe the initial wake; the user may describe the line geo- 
metry in as much detail as he likes (1). The program inter- 
rogates each line in relation to a system of vertical WAKE-GRID 
PLANES, Figure 11. The user specifies the location of these 
planes in accordance with the anticipated track of the wake; 
because the wake lines are represented by straight line segments 
between the wake-grid planes, regions where the wake lines are 
highly curved require an increased density of wake-grid planes. 

Wake panels are contained between adjacent wake-grid planes 
and adjacent wake lines. The wake panel geometry is evaluated in 
the same way as the surface panels (Subsection 3.1.2) except that 
the mean plane is constrained to lie in the upstream edge (i.e.. 
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plane containing panel normal vector , 

CONTROL POINT AND MID - POINT OF SIDE 
AT PATCH EDGE 



Figure 10. Panel Neighbors Across Patch Edges. 



TYPICAL WAKE-6RID PLANES 
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Figure 11. Wake-Grid-Plane Scheme. 








Figure 12. Wake Arrangement 





side 4) and to contain the mid-point of the downstream edge (side 
2). In that way the projected wake panels are cleanly attached 
to the surface separation line even when the panels are skewed. 

The wake geometry supplied by the user is used only for the 
initial solution. The program includes an iteration cycle in 
which the straight wake-line segments (i.e., panel edges 1 and 3) 
are aligned with the computed flow direction (see Subsection 
3.5) . 


3.2 Matrix of Influence Coefficients 

Equation (10) in Section 2.0 is a Fredholm integral equation 
of the second kind for the unknown velocity potential distribu- 
tion over the body surface. A velocity potential distribution 
which satisfies Eq., (10) at every point of the body surface 
provides an exact solution for the inviscid, incompressible po- 
tential flow about the configuration surface through Eq. (13). 
In order to obtain a numerical solution the equation is satisfied 
only at a finite number of points on the surface; i.e. f at one 
point (the control point) on each panel. Also, the doublet and 
source singularity distributions are assumed to be constant on 
each panel and the surface integrals in Eq. (10) are performed in 
closed form over each panel. These integrals, therefore, provide 
the influence coefficients per unit singularity strength for each 
panel. The total surface integrals then become summations over 
all panels and Eq. (10) is transformed into a set of simultaneous 
linear algebraic equations for the doublet strength on each 
panel . 

N S N w N s 

IX C JK> + Z V C JL> + Z °K B JK ' ° ! J = 1 ' Ns 

K=1 L=1 U K=1 


( 20 ) 

where y and o are given by Eqs. (11) and (12), respectively, and 
C JJ = " 2it * 

N s and N w are the number of surface and wake panels, respec- 
tively. The quantities, B JK and C JK , are the perturbation 
velocity potential influence coefficients for the constant source 
and doublet distributions, respectively, on panel K acting on the 
control point of panel J; i.e., from Eq. (10). 
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V (-) 
r 


dS 


( 21 ) 


where n = n K is a constant for a flat panel, K, 
length of the vector from the surface elements, dS, 
the control point of panel J. 
to the surface point. 


and r is the 
of panel K to 
The gradient is taken with respect 


Equation (20) is applied to all panels, J, that occur on 
"Dirichlet" patches (IDENT=1 or 2); i.e., on patches that are 
part of a "thick" boundary enclosing an inner volume. If the 
control point, J, is on a panel on a "Neumann" patch (IDENT=3), 
i.e., part of an open surface, then a discretized form of Eq. 
(19) is used. 


Z % e jk> + Z % e jl> 

K=1 L=1 


+ £ °K D JK + "j 

K=1 


V 



0 


( 22 ) 


where 


D 


JK 


n 


V 

a 


JK 


and 



(23) 


are, respectively, the normal velocity influence coefficients for 
the source and doublet distributions on panel K acting at the 
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control point of panel J. The velocity vector influence coeffi- 
cients for the source and doublet are, respectively. 


V 

a 


JK 


- // V7 ,ds 

Panel K 


and 



- // V"K 

Panel K 


V(^) )dS 


(24) 


with the gradients, V jr being taken with respect to the coordi- 
nates of point J. 


Note that in generating Eq. (19) for Neumann patches, the 
source terms disappeared with the assumption of no jump in normal 
velocity through the surface. However, for the general problem 
involving a combination of thin and thick components the source 
term has been left in Eq. (22) to supply the contribution from 
source panels on thick boundaries. 

The influence coefficients for the velocity potential and 
velocity are developed separately below. In each case the coef- 
ficient is determined completely by the geometric details of 
panel K and the coordinates of the control point. 


3.2.1 Velocity Potential Influence Coefficient 

3. 2. 1.1 Doublet Contribution 

Consider the doublet velocity potential influence coeffi- 
cient from Eq. (21). After evaluating the gradient, this becomes 



Panel K 


(25) 
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The integral is evaluated by considering each side of the 
(projected) panel in turn, each side having a pair of semi- 
infinite strips of opposing singularity strength equal to half 
the panel value, Figure 13. The strips lie in the mean plane of 
the panel and are aligned with the panel's & axis. With the 
convention adopted for the panel sides (i.e., the side order is 
anticl ockwise when viewed from the flow field) the singularity 
strength is given a positive sign to the left of the side and a 
negative sign to the right when looking in the direction of the 
side. Figure 13. When all four sides of the panel (or for any 
closed polygon, in fact) have been treated in this way then the 
contributions from all the overlapping strips cancel outside the 
panel and reinforce for unit singularity value over the panel. 

Consider the contribution from one side. The kernel can be 
simplified by combining corresponding pairs of positive and nega- 
tive elements, d£, dy, located at equal distance, E , , from a point 
(x,y) on the side (Figure 14); i.e.. 


C 


JK 


Y b - 

•'// 

x 7— r = n 




<H dy 


(26) 


where 

r = R(n) + ££ 

1 

r = R(n) - C*- 

2 

R(n) = a - ns 

n = (y - y a ) / <Yb ~ y a > 


where y a , y^ are the beginning and end of the projection of the 
side on the m axis and a is the position vector of the J th 
control point relative to the start of this side (Figure 14). 

The final expression for the doublet velocity potential 
contribution from side i of the panel is (using variables which 
closely follow the computer coding): 


- 1 

C JK =tan (RNUM/DNOM } 

l 


(27) 
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Figure 14. Evaluation of the Surface Integrals for One Side of a Panel. 




where 


RNUM = SM * PN * (B * PA - A * PB) 
DNOM = PA * PB + PN 2 * A * B * SM 2 


PN P JK * n K 

P JK = Rc J “ Rc k 

A = lal 

B = I b I 

PA = PN 2 * SL + A1 * AM (= a • (t x (a x s )}) 

PB = PN 2 * SL + A1 * BM = PA - A1 * SM 

SM = s • m 

SL = s • l 

AM = a • m 

BM — b • on 

Al = AM * SL - AL * SM (= n • (S x a)) 

b is the position vector of the J th control point 

relative to the end of this side. Figure 14 (i.e., b = 
a - s) 


The subscript, i, has been omitted from the above variables; a, b 
and s are constructed for each side in turn, noting that a. = 

b i-l* 

There is a limiting condition as PN — the projected height of 
control point J from the mean plane of panel K — goes to zero: 

= ± tt within the strip (DNOM<0) 

= ±\ on the edge of the strip (DNOM=0) 

= 0 outside the strip (DNOM>0) 

The positive sign occurs for points to the right of the side 
and negative to the left. If PN •+■ 0", all the above signs on tt 
are reversed. 

The total velocity potential influence coefficient for the 
unit doublet on the panel is the sum for all four sides, from Eq. 
(27) . 


Lim (C- rv ) . 

+ K i ' 

PN^O 


37 


(28) 


C 


JK 


E 

■; = 1 



The four arc-tangent functions in this summation can be 
combined into a single arc tangent; however, extreme care is 
required to ensure the resultant angle is in the correct quadrant 
in the range, ± 2 tt . 

For plane of symmetry and/or ground effect problems the 
contribution to Cjk from image panels (having the same singulari- 
ty strength as the real panel) is evaluated simply by repeating 
the calculation for the real panel influencing the images of 
control point J. This is possible because the geometrical rela- 
tionship between the image of panel K and the point J is identi- 
cal to that of the actual panel K and image of point J (Figure 
15). Hence, for a symmetrical problem in ground effect there 
would be four contributions for each panel. Figure 15. This 
treatment is more convenient than considering the image panels 
directly. The vertical plane of symmetry is the y = 0 plane in 
the G.C.S., while the ground plane is the z = 0 plane. Thus, for 
ground effect problems the configuration height and orientation 
above the ground require care during input (1). 


The above expression for Cj K is used when panel J is close 
to panel K. When panel J is remote from panel K, i.e., in the 
far field, then a far-field formula is used under the assumption 


n • it 

that the quantity, — , in Eq. (25) is essentially constant 

r 

over the panel. Then r = P JK and 


C 


JK 


PN t „ * AREA 

U 1\ J\ 



(29) 


where 


PN JK = 
n K = 

area k = 


P JK ' n K 

the unit normal to panel K 

the area of panel K (see Subsection 3.1.2) 

lPj K l, is the distance between the control points, 
and J 


K 


Equation (29) is the contribution from a point doublet at 
Rc k with orientation n K and is used when P™ > RFF * (SMP + 
SMq) k . The quantity, (SMP + SMQ)r/ is the "size" of panel K (see 
Subsection 3.1.2) and RFF is the user-defined, far-field radius 
factor which has a default value of 5 in the program. This value 
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xc. 






was established following initial test cases on a wing and on a 
wing-body configuration. Use of the expression gives a signifi- 
cant reduction in computing time; however, a user may use the 
near-field equation (Eq. (25)) throughout the configuration if he 
wishes simply by giving a very large value for RFF in the input 
file (1). 


3. 2. 1.2 Source Co ntribution 

The velocity potential influence coefficient for a unit 
uniform source distribution on the panel, K, is, from Eq. (2 , 



Panel K 


The integral is treated for each side of the panel as in the case 
of the doublet. Thus, for one side 


B 


JK i 


= k 


Yb 

I 


!j_ i I 


i r i 


r 2 i 


d£, dy 




where r , r , y , y are defined under Eq. (26). This becomes: 

1 2 ci b 


B = A1 * GL - PN * C JK _ 


(30) 


where 


and 


GL 

S 


I LOG 
s 

I s I 


A + B + s 
A + B - s 


Al, PN, C jKi/ A and B are defined in Eq. 
second LOG term arises in the above integra 
tribution from all sides of a closed polygon 


(27). (Note that a 
1 but its total con- 
cancels exactly.) 
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Again, the total source velocity potential influence coeffi- 
cient is obtained by summing the contribution from all four sides 
of the panel; i.e.. 


4 



i=l 


(31) 


Symmetry and/or ground-effect treatment is the same as that 
described for the doublet contribution. 

The expression for B JK for points that are in the far field 
from panel K is, using a similar approach to that for the doublet 
(Eq. (29)), 


area r 

Bjk ' 


(32) 


With the above expressions for C JK , B JK substituted in Eq. 
(20), the zero internal perturbation formulation is similar to 
that given by Morino (19), who applied Green's Theorem directly 
to the external flow. 


3.2.2 yelocitv_Vector Influence Coefficients 
3. 2. 2.1 Doublet Contribution 

The velocity induced at control point, J, by a constant 
doublet distribution of unit value on panel K is, from Eq. (24), 


if ^ ’ V< ^ >,dS 

Panel K 


(33) 


The kernel is one of the terms resulting from the expansion 
of a triple vector product. Hence, V can be written 

y JK 


JK 


-// 


{(n x V ) x 7(1) + nV., * V(-)}dS 

J' r J r 


But, V a • V(±> 


Panel K 

= 0 . 
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Hence , 


JK 



(n x Vj) x V(i)dS 


(34) 


Panel K 


or 


JK 




V (— ) x dr by Stokes' Theorem, 

r 


Thus, 


JK 


-p-i? 


dr 


(35) 


The integral can be performed along the straight line of 
each side of the panels Figure 16. Hence, the contribution from 
side i is 


a x b 


y a x b • a x b 

JK • 

1 


* (A + B) 


/A * B 
\ A 


- a • b 


B 


This expression is singular when the point, P, is in line 
with the side (a x b = 0). However, 


a x b • a x b = (A * B - a • b) (A * B + a • b) . 


And so an alternative expression is 


a x b * (A + B) 

u - A*B* (A*B + a*b) 

m jk. 

1 


(36) 


This becomes singular only when the point, P, lies on the 
side. The singularity is avoided by applying a finite-core model 
on the vortex line. (The core radius can be set by the user 
( 1 ) .) 

The total velocity coefficient induced by the unit doublet 
on the panel is the sum of Eq. (36) over all four sides. 
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(37) 
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JK 


£ 

-i - 1 


JK 
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For points in the far-field, the integrand in Eq. (33) is 
treated directly as a point doublet; i.e. f 


1 = - V T (n__ • V(^)) AREA 

UjK J K r K 


or 


J = - AREA * V 

y jK K J 


(^) 


Thus, 


JK 


AREA k * (3 * PN jk * P JK - P JK * «K 


}/P 


JK 


(38) 


The above variables are defined in Eq. (29) . 

Cases with a vertical plane of symmetry and/or ground effect 
are treated as before by evaluating the influence of the 
panel at each i ma g e point rather than vice versa. With this 
approach the velocity vector contribution must have its component 
normal to the symmetry plane reversed in sign in order to obtain 
the actual contribution from the image panel acting on the real 
control point. 

The normal velocity influence coefficient, E JR , in Eq. (23) 
can now be evaluated; i.e.. 


E 


JK 


n. 


JK 


Where the velocity influence, V , is taken from either Eq. (37) 

m jk 

or Eq. (38) depending on the distance of the point, J, from panel 

K. 


3. 2. 2. 2 Source Contribution 

The velocity induced at control point, J, by a constant 
source distribution on panel K is, from Eq. (24), 


JK 


-// 


(^)dS 


Panel K 
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and taking the gradient with respect to the coordinates of point 
J, this becomes 



Panel K 


(39) 


The integral is treated for each side of the panel as in the 
case for the velocity potential evaluation. From side i f the 
contribution is 

Vo = GL * (SM * A - SL * m) + C n (40) 

JK . 

l l 

Where GL is defined under Eq. (30), and SM, SL and C JK . are 
defined under Eq. (27). The total for the panel is the sum" over 
all four sides. 

4 

V0 JK ’ T, V °JK. <41) 

i=l 1 


For points in the far field, the panel's source distribu- 
tion is concentrated in a single source at the panel center. In 
this case the source velocity coefficient becomes (from Eq. (39)) 


Vo 


JK 


area k * P JK 


JK 


The above variables are defined in Eq. (29). 


(42) 


Cases with a vertical plane of symmetry and/or ground effect 
are treated as before by evaluating the influence of the real 
panel at each imas£ point rather than vice versa. With this 
approach the velocity vector contribution must have its component 
normal to the symmetry plane reversed in sign in order to obtain 
the actual contribution from the image panel acting on the real 
control point. 


Finally, the normal velocity influence coefficient for the 
source term in Eq. (23) can be evaluated; i.e.. 


The source veloc 


Eq. (41) or Eq. (42), 
panel K. 


D — n • ¥o 

JK J JK 

ity coefficient, V a , is taken from either 

JK 

depending on the distance of point J from 
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3.2.3 Matrix Assembly 


The matrix equations, Eq. (20) and (22), can now be formed 
by assembling the singularity influence coefficients for all the 
surface and wake panels in the configuration acting at the con- 
trol point of each surface panel in turn. In the VSAERO program 
the surface and wake contributions are computed separately: the 
surface matrix of influence coefficients, once computed, remains 
fixed throughout the treatment of a configuration while the wake 
contribution varies with each wake-shape iteration. The wake 
contribution is, in fact, recomputed and combined with the wake- 
shedding panel influence coefficients after each wake-shape cal- 
culation. 

Consider the wake influence in Eq. (20); i.e.. 


W 

yi u w c jl 

L=1 


Now from the boundary condition imposed on the wake (Eq. 
(7)) the value for y w is constant down each wake column and has 
the value of the jump in total potential at the start of the 
column (from Eq. (10)); i.e., for the Mth column, 


U W„ “ °KWPU X , ” ^KWPL.. 
M M M 


is constant down the column; KWPU^, KWPL^ are the upper and 

lower surface panels shedding the wake column (Subsection 

3.1.4). 


Thus , 


y w„ " y KWPU M 
M M 


y KWPL.. 

M 


+ <j> 


°°kwpu m 

M 


°°kwpl m 

M 
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And the wake influence contribution can be written 


N w 

NWCOL 

y w c jl = 

£ 

L=1 

M=1 


I 


PHC JM * (y KWPU y KWPL M ) 

M 


+ PHC 


JM 




KWPU 


- <f> 

' n 


M 


KWPL, 


( 44 ) 


M 


where 


NWP 


M 


PHC 


JM 


L=1 


'JL 


NWCOL = Total number of wake columns 


and 


NWP.. = 


M 


Number of wake panels on the M th column 

cient^*^ the s 4 mmation of the doublet influence coeffi- 
subtract<Td fro?!? 11 ' 1,e Y PI l c JMf is therefore added to and 
Danpi inn d from * he associat ed upper and lower wake shedding 
P n luences (KWPU^, KWPIj^), respectively. 

, i® tempting to discard the last term of Eq. 

the quantity ( ^ - <J> ro ) vanishes in the 

kwpu m kwpl m 

Drnarh ^ We . r 1 contro1 points on the wake-shedding panels ap- 

P^actical oanel d^nJ^ 6 in # .S^ gh densit y Paneling. However, with 
? n v^. 31 panel densities, this quantity will rarely be zero and, 

ct, Youngren et al. (20) noticed that a serious error in 
circulation results if the term is omitted for thick trailinq 

The term 1S ' in fact ' readily evaluated and is combined 
with the source terms to form the known "right-hand side" of the 
equations. 


(44) because 
limit as the 
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3.3 r i x Solver 

With the matrix of influence coefficients formed, the solu- 
tion routine can proceed. VSAERO employs either a direct or an 
iterative matrix solver depending on the number of unknowns 
(i.e., number of panels). The direct solver is Purcell's vector 
method (21) and is used when the number of panels is less than 
320. When the number of unknowns is greater than 320, then a 
"blocked controlled successive under- or over relaxation \IA) 
technique is used based on the Gauss-Seidel iterative technique. 

The iterative technique proceeds by forming inverses of 
smaller submatrices centered around the diagonal and by mul i 
plying the residual error vector with the stored inverse of each 
submatrix. The submatrix blocks are formed automatically within 
the code on the basis of complete columns of panels. An option 
is provided for the user to specify his own blocks when the 
arrangement of patches and panels becomes difficult in awkward 
geometrical situations. Placing two strongly interacting panels 
in an order which prevents their inclusion in the same oc oan 
lead to poor convergence and sometimes divergence in the blocked 
Gauss-Seidel routine. 


3.4 Qn-Bodv Analy sis 

The matrix solver provides the doublet (i.e., the perturba- 
tion velocity potential) value for each surface panel in the 
configuration. These values and the source values are now as- 
sociated with the panel center points in the evaluation of the 
velocity on each panel. 

The local velocity is: 


V = - 7$ 


The perturbation velocity, v, is evaluated in the panel's 
local coordinate system; i.e., 

v = VLi, + VMm + VNn (46) 

The normal component is obtained directly from the source 
form, Eq. (12); i.e., 

VN = 4tto (47) 

The tangential components, VL, VM, are evaluated from the 
qradient of y, assuming a quadratic doublet distribution over 
three panels in two directions in turn, Figure 17. Thus, when 
evaluating the velocity on, say, panel K, the four immediate 
neighboring panels, Nl, N2, N3, N4, are assembled together with 
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ure 17. Evaluate 
Availab] 


the neighboring sides of those panels, NS1, NS2, NS3, NS4. This 
information is obtained directly from stored neighbor information 
( (NABOR(I ,K) , NABSID (I ,K) , 1 = 1,4), K = 1,NS), see Subsection 3.1.3. 
Surface distances are constructed using the SMP, SMQ information 
for each panel (Subsection 3.1.3), noting that the distance from 
the center of panel K to the midpoint of either side 1 or side 3 
is SMQr and, similarly, the distance to the midpoint of either 
side 2 or side 4 is SMPk* The additional distance from the 
midpoint of a side of panel K to the midpoint of the neighboring 
panel of that side depends on the NABSID value. For example, the 
additional distance to the neighbor, Nl, from side 1 of panel K 
is SMQ Nl if NABSID ( 1 , K ) is 1 or 3 and is SMP Nl if NABSID(1,K) is 
2 or 4; the former is always the case when the neighbor is in the 
same patch as panel K (Figure 17) and the latter will occur when 
the neighbor is on a different patch with a different orienta- 
tion . 


Consider the regular case where the neighboring panels are 
available and are on the same patch as panel K. Then the deriva- 
tive of the quadratic fit to the doublet distribution in the SMQ 
direction on panel K can be written 


where 


DELQ = (DA * SB - DB * SA)/(SB - SA) 


SA = 

- ( smq k 

+ SMQ, 

SB = 

SMQ k + 

SMQ 

N3 

DA = 

<l, Nl " 

u k >/sa 

DB = 

<1J N3 ' 

u k >/sb 


(48) 


A similar expression based on SMP rather than SMQ and using 
neighbors, N4 , N2, instead of Nl, N3, gives the gradient DELP in 
the SMP direction. 

DELQ and DELP provide the total gradient of y and, hence, by 
resolving into the l and m directions, the perturbation tangen- 
tial velocity components (Eq. (46)) can be written 

VL = - 4 it * (SMP * DELP - TM * DELQ) /TL 


and 


VM = — 4 7T * DELQ (49) 

The vector, TL, TM, with modulus, SMP, is the center of side 
2 relative to the panel center expressed in the local coordinate 
system of panel K. 
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f vl he c u omple J te set of Panel neighbors is not always available 

c?Li : i 1 e P H ab HJ e 1 .K derivative °P erat ion. Neighbor relationships are 
cancelled by the program across each wake shedding line to pre- 

vent a gradient being evaluated across a potential jump* The 
program also terminates neighbor relationships between panels 
when one has a normal velocity specified and the other has the 
default zero normal velocity. in addition, the user can prevent 
neighbor relationships being formed across a patch edge simply by 
setting different ASSEMBLY numbers for the two patches and, 
finally, the user can actually terminate neighbor relationships 
by-panel basis in the input. This is applied wherever 
tnere is a potential jump (e.g., where a wake surface butts up to 
the fuselage side) or where there is a shape discontinuity in the 
flow properties (e.g., a sharp edge) or where there is a larqe 
mismatch in panel density from one patch to another. 

When a panel neighbor is unavailable on one or more sides of 
a 4 T P ^u e1 ' . ? n the dou blet gradient routine locates the neighbor 
of the neighbor on the opposite side to complete a set of three. 
For example, if neighbor N2 is not available. Figure 18, the 
neig or on the opposite side is N4 with adjacent side NS4. 
Going to the far side of panel N4, i.e., away from panel K, the 
neighboring panel is NABOR (MS4 ,N4) with adjacent side 
NABSID (MS4,N4) , where J e 


MS4 = NS4 ± 2; 


1 « MS 4 ^ 4 


^ this pane ^ available then the set of three is com- 
' A, gradient now being taken at the beginning rather than 
the middle of the quadratic. If the panel is not available then 
the program uses simple differencing based on two doublet values, 
lit no panel neighbors are available, then the program cannot 
evaluate the tangential perturbation velocity for the panel K.) 

With the local velocity known, the pressure coefficient is 
evaluated on panel K: 


C 


PK ■ 


1 


v 2 „/v 2 

K <» 


(50) 


The force coefficient for a configuration or for parts of a 
configuration are obtained by summing panel contributions. The 
contribution from panel K is (Figure 19) 


AF k 

— = - Cp K * ARBA k * n K + C f * AREA k * t k (51) 

K 


where 

^oo * s the reference dynamic pressure 

C P K ts the pressure coefficient evaluated at the panel 

center, Eq. (50) 
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Figure 19. Force Contribution from Panel 



AREA„ is the panel area projected onto the panel's mean 
plane 


n is the unit vector normal to the panel's mean plane 

K 

t is the unit vector parallel to the panel's mean plane 

K and in the direction of the local (external) 
velocity; i.e., t R = V k /I\M 

C is the local skin friction coefficient 

f K 


The force vector for a patch is then 



(52) 


where Ki, are the first and last panel subscripts in the 

patch. The force coefficient vector is 


K 0 

2 AFr 

C F * V q “ Sref = q °° /Sref 

P P K= Kl 


where S is the reference area specified by the user. 

REF 

The components of the force coefficient vector are printed 
out in both wind axes, C L , C D , 0$ , and body axes, C x , C y , C z , 
(i.e., the global coordinate system). Figure 20. 

For symmetrical configurations the total force vector is 
doubled prior to printout. This means that the printed side 
force value, i.e., the component normal to the plane of symmetry, 
is actually twice the side force on one side of the configuration 
rather than being cancelled to zero for the total configuration. 


The skin friction coefficient, Cf , in Eq. (51) is zero for 

K 

the first time through the potential flow calculation. During 
the viscous/potential flow iteration cycle the Cf values are 
computed along the surface streamlines passing through user- 
specified points. The program then sets the local Cf value on 
each panel crossed by a streamline. Panels not crossed by 
streamlines have C f values set by interpolation through local 
panels that have been set, provided these are in the same patch. 
Patches having no streamlines crossing them are left with zero Cf 
values on their panels; consequently, it is up to the user to 
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' Global Coordinate 
System (Body Axes) 


Figure 20 . Global and Wind Axis Systems. 
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adequately cover the regions of interest with streamline requests 
(1). If a patch is folded (i.e., a wing) and/or there is a 
dividing stream surface approaching the patch, then streamlines 
should be requested on both sides of the attachment line. 


The contribution to the 
panel is: 


loment coefficient vector from the 


XAF K 


where R K is the position vector 
panel, K, relative to the moment r 
reference semispan for the x and z 
and yawing moments); and L=CBAR (t 
component of the moment (pitching 
specified by the user (1). 


of the control point on the 
eference point; L=SSPAN, the 
components of moment (rolling 
he reference chord) for the y 
moment). CBAR and SSPAN are 


The summation is ca 
ponents and assemblies, 
acting on user-specified 


rried out on the basis of patches, com- 
Summation of the forces and moments 
sets of panels can also be obtained. 


Patches with IDENT=1 have an additional output of local 
force and moment contribution from the summation down each patch. 
This is useful for "folded" patches representing wing-type sur- 
faces. 


Patches with IDENT=3 include both upper and lower surface 
velocity and pressure output. On these patches the doublet 
gradient (Eq. (49)) provides the jump in tangential velocity 
component between the upper and lower surfaces. The program 
first computes the velocity (i*e., excluding the loca 

panel's contribution) at the panel's center then adds to this one 
half of the tangential velocity jump for the upper surface and 
similarly subtracts for the lower surface. 
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3.5 Wake-Shape Calc ulation 


With the panel singularity values known, the velocity vector 
can be computed at any point in the flow field using the panel 
velocity influence coefficients and summing over all surface and 
wake panels (Eq. (18)). In the wake-shape calculation procedure 
the velocity vector is computed at each wake panel corner point 
and each "streamwise" edge; i.e., sides 1 and 3 of a panel. 
Figure 11, is aligned with the average of its two end velocities. 

The wake shape calculation is performed in a marching pro- 
cedure, starting at an upstream station and progressing through 
the set of vertical wake-grid planes. Figure 11. In each wake- 
grid plane the full set of (streamwise) wake lines intersecting 
that plane is first assembled and the full set of velocities 
computed at that station. In assembling the set of wake lines 
for treatment, certain lines can be left out, e.g., lines on 
wakes identified as rigid by the user (1); individual lines that 
have been identified as rigid, for example, the line where the 
wing wake butts up to the fuselage side; and lines that are 
coincident with other lines that are already in the set. 

The step from one wake-grid plane to the next is performed 
for all lines before another set of velocities is computed. 
Thus, for a given set of points, RW ± , M ; i = l, NL M , in the M th 

wake-grid plane, the set of points in the next wake-grid plane 
is: 

j J j j 

RW i,M+l “ ""i.M + 5 i * DX/IV x .l; i-l, NL m (55) 

where NL^ is thenumber of wake lines being moved in the M th 
wake-grid plane; V? is the mean velocity across the interval for 
the i z line in the j th predictor/corrector cycle. 

_ 3 i J 3-1 \ 

V i - ( V i,M + V i,M + J' 2 (56) 


J is the total number of predictor/corrector cycles in use (de- 
fault is 2). 

o _ j j-1 

. ^ 1^+1 is set of V i /M ' and v i,m+ 1 is coin P uted at the point, 
RW i,M+l • 

In each predictor/corrector cycle the complete set of points 
in each wake-grid plane are moved and the current movement vector 


d 


j j-1 

RW, - RW „ 

i,M i,M 
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is applied to all the points downstream on each line, before the 
next set of velocities is computed. A limiting movement can be 
applied as a damper in difficult situations. 

When all the wake-grid planes have been processed, the wake 
panel parameters are reformed using the new corner points. The 
modified set of influence coefficients can then be computed for 
the next solution. 


3.6 Boundary Layer Analysis 

The solution from the on-body analysis (Subsection 3.4) can 
be transferred to the boundary layer analysis. Two integral 
boundary layer analyses are provided in the VSAERO program: the 

first method, which is described in Reference 23, provides the 
boundary layer characteristics along computed (external) surface 
streamlines. The location of a point (panel) on each streamline 
is provided by the user who is, therefore, responsible for en- 
suring that a family of streamlines adequately covers the regions 
of interest. The method, which is applicable for wings and 
bodies, includes the effects of surface curvature and streamline 
convergence/divergence under the assumption of local axisymmetric 
flow. The surface streamline tracking procedure is described in 
(8) and (24) and is based on a second-order surface stream 
function formulation. The streamline boundary layer calculation 
can be expected to break down in regions of large cross flow. 
However, the procedure has given surprisingly good guidance for 
the onset of separation for a wide range of practical problems. 

The second method, which is installed in the VSAERO program 
and which is described in Reference 2, includes a cross-flow 
model. This method is an infinite swept wing code and is applied 
in a strip-by-strip basis across wing surfaces, with the user 
specifying which strips are to be analyzed. 

The boundary layer module is fully coupled with the VSAERO 
inviscid code in an outer viscous/potential iteration loop. 
Figure 3. In further iterations, the displacement effect of the 
boundary layer is modelled in the potential flow calculation 
using the source transpiration technique; i.e., the boundary 
layer displacement term appears in the external Neumann boundary 
condition equation as a local non-zero normal velocity component 
(Eq. (12)), and is thereby related directly to source singulari- 
ties in this formulation. The boundary layer calculation also 
provides a separation point on each streamline, thereby defining 
a boundary of separation for a family of streamlines. The cri- 
terion for separation is currently a vanishingly small skin 
friction coefficient (in the direction of the external stream- 
line) . 
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3,7 Qff-Body Analyses 


hnnn^ en i^ e iterative loops for the effects of wake shape and 

for field v Ye i r ar ^ e complete ' the VSAERO program offers options 
tor field velocity surveys and off-body streamline trackina 

off-bodv^eln'r " hlch ara discussed separately below, both use the 
off body velocity calculation, Eq. (18), based on a summation of 

onset °flow? Utl ° nS fr ° m a11 surface and wake Panels added to the 


3*7.1 Veloci ty Survey 

body velocity calculations are performed at user selec- 

scan^l ines* Ihl ™* th ® points are assembled along straight 
f J n ! " es ‘. These lines in turn may be assembled in planes and 

volun,es * The shape of the scan volume is deter- 
able m Parameter, MOLD. Two options are currently avail- 


MOLD 1 Allows a single point, points along a straight 
line, straight lines within a parallelogram 
within a parallelepiped. 


or 


MOLD=2 


Allows points along radial lines in a cylindrical 
volume. 


A number of scan boxes can be specified, one after the 
other. A MOLD=0 value terminates the off-body velocity scan. 

4 . Th , e , \ oc . ation of scan lines within the scan boxes can be 

generated^ lo th f, inpU *T The default option is equal spacing 

ff?* J e e d al ?" g the sides of the box - The location of points 

the dPfafn? 11 1 l™ Can be controlled also in the input, but P again, 
tne default option is equal spacing. 

an line * n . tersects the surface of the configuration, 

. . tlon routine locates the points of entry and exit 

each U ?nVor e sai \ face Paneling. The scan line points nearest to 
each intersection are then moved to coincide with the surface 
Points falling within the configuration volume are identified by 
e routine to avoid unnecessary velocity computations These 

controls \1ip 1 ra 9 i 9 i ed ^ n ^ he printout * The input parameter, MEET, 
ntrols the calls to the intersection routine: MEET=0 (default) 

makes the routine active, while MEET=1 switches it off; clearly, 
if the selected velocity scan volume is known to be out s iH P ^ 

co°m n L\ 9Ur w a H tl0n r 1Ume ' t hen th6re is no need for the program to 
compute where the scan lines intersect the configuration surface. 

The velocity routine, VEL, computes the velocity vector at 

wlke ^a'neSs 7 th 1 CO /^ r ibutions from all surface panels, 

wake panels, image panels (if present), onset flow, etc. The 
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routine includes near-field procedures for p oints 

i-hat are close to the surface panels or wake panels. The near 
field routine involves a lot of extra computation; if, therefore, 
is known that a scan box is well clear of the configuration 
surface and wake then it is worthwhile turning off the near field 
routine to avoid unnecessary computation. This 1S J ontro11 ® .X 
the input parameter, NEAR; NEAR=0 (default) keeps the near field 

routine active while NEAR=1 switches it off. 


3.7.2 of f-Bodv Streamline Pro cedU££ 

The numerical procedure for calculating streamline paths is 
based on finite intervals, with a mean velocity vector being 
calculated in the middle of each interval. Figure 21. The 
velocity calculation point, RP, is obtained by extr ^P°^ tlon t g 
the two previous intervals on the basis of constant rate of 
change in the velocity vector direction; i.e., 

RP = R n + * 5 S n t (5?) 

where R is the previous point calculated on the streamline, s 
Is the present interval length, and t is a projected tangent 
vector . 


where t = , etc. V . is the velocity calculated in the 

wnere c r _ 1 |V n-1 1 n-1 

middle of the previous interval. 

RP does not necessarily lie on the stream line path. 
calculated velocity, V n , at RP, is used to evaluate the 


point on the streamline; i.e. f 


R n+1 = R n + 


s V 

n n 


An automatic procedure is included in the routine which 
changes the interval length, s n , in accordancewiththechangein 
tanoent direction. If a large change in direction is calculated, 
the 9 value of s n is decreased and the calculation is repeated 
until the change in direction is within a specified amount. If 
the change in direction is smaller than a certain amount, then 
the interval length for the next step is increased. 
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STREAMLINE PATH 


/ 



VELOCITY CALCULATED 
IN MIDDLE OF EACH INTERVAL 
GIVES LOCAL TANGENT VECTOR 


Figure 21 . 


Streamline Calculation. 


4.0 DISCUSSION OF RESULTS 


Correlation studies discussed below include exact and nomi- 
nally exact solutions of a number of basic test cases, some of 
which had proved difficult for earlier low-order panel methods. 
In addition, some comparisons are made between calculated and 
experimental data. 


4.1 Swept Wing 

Figure 22 shows the chordwise pressure distribution at u = 
0.549 on a thin, swept wing at 5° incidence. The wing has a mid- 
chord sweep of 30°, aspect ratio 6, taper ratio 1/3, and a NACA 
0002 section. A nominally exact datum solution (25) for the case 
was computed by Roberts' third-order method using a 39 
(chordwise) by 13 (spanwise) set of panels. The present calcula- 
tions, based on the <*>^ = 0 formulation, use 10 spanwise panels 
with both 80 and 20 chordwise panels on a "cosine" spacing, 
giving increased panel density towards leading and trailing 
edges. Roberts' datum paneling is heavily weighted towards the 
leading edge: he has, in fact, 4 more panels ahead of the most 
forward panels in the present 80 cosine spacing case. 

The comparison in Figure 22 indicates that the low- and 
high-order methods give essentially the same solution where the 
control point density is the same. The present calculations 
using 20 panels chordwise deviate from the datum solution near 
the leading and trailing edges where the control point densities 
are different. Good agreement is restored in those regions by 
the 80-panel case which, at the trailing edge at least, has a 
similar panel density to the datum case. (Note: not all the 80- 
panel pressure values are plotted in the central region of the 
chord in Figure 22.) 


4.2 Wing with Strake 

The fact that surface velocities are obtained by numerical 
differentiation in the present method might lead to inaccuracies 
in awkward situations. So far, no serious trouble has been 
experienced. Any problem that might arise with the derivative 
scheme (e.g., in a region of mismatched paneling) can usually be 
alleviated by cutting panel neighbor relationships at that sta- 
tion , thus forcing a forward or backward differencing (Subsec- 
tion 3.4). One situation that could be difficult and which has 
caused a problem for other methods is the evaluation of the 
spanwise velocity component, Vy, in the neighborhood of the kink 
on a wing with strake. In the case considered in Reference 25, 
the leading-edge kink occurs at 25% of the semispan and the 
leading-edge sweep is 75° inboard and 35° outboard. The trail- 
ing-edge sweep is 9°. Again, the wing section is a NACA 0002 and 
incidence 5 °. 




22. Comparison of Calculated Chordwise Pressure Distri 
butions on a Thin Swept Wing. NACA 0002 Section; 
a = 5" 


The chordwise distribution of Vy is plotted in Figure 23 at 
stations n = -219 and .280; i.e., just inboard and just outboard 
of the kink, respectively. 

Two hi ghe r - o r de r datum solutions are plotted (25) from 
Roberts' third-order method and from Johnson and Rubbert's (17) 
second— order method. These datum solutions agree on the outboard 
distribution but disagree by a small amount on the inboard cut. 
The disagreement might be caused by the different ways the inter- 
polation is set up across the kink. 

The present calculations are in close agreement with the 
datum solutions and favor the second-order solution near the 
upstream part of the inboard cut. The case was run as part of a 
"shakedown" exercise of the code and used 946 panels in a 4 
(chordwise) by 20 (spanwise) array with 66 panels on the tip— edge 
surface. The chordwise paneling was on a "sine" distribution 
with density increasing towards the leading edge. Such a high 
number of panels was probably not needed for this case, but the 
calculation took approximately 2^ minutes on a CDC 7600. 


4.3 Nacelles 

In the past, low-order solutions for internal flows have 
suffered badly from "leakage" problems in between control points. 
The present code was, therefore, applied to two nacelle cases. 
The first of these is a simple through-flow nacelle formed by 
rotating a NACA 0010 about the nacelle axis such that the minimum 
area is one fourth of the inlet and exit areas. Figure 24. This 
example, taken from Reference 18, shows how the earlier, first 
generation, low-order analysis, because of the leakage, cannot 
match the internal flow calculated by the higher-order methods. 
However, the low-order, VSAERO solution (with the internal poten- 
tial set to zero) is in close agreement with the high-order 
solutions and evidently does not have this leakage problem. 

The second nacelle case considered is also an open nacelle, 
but with an exit diameter 0.3 of the chord. The nacelle surface 
is generated by a NACA 0005 section tilted leading edge out by 5° 
from the x-axis. 

Two higher-order datum solutions (25) are shown in Figure 
25. for the chordwise pressure distribution. These are due to 
Hess and also Roberts. 

Roberts' solution shows a slightly lower pressure inside the 
duct compared with Hess's solution. The present calculations, 
based on the = 4^ formulation, use 10 panel intervals around 
half the circumference and with two chordwise distributions, 26 
cosine and 48 sine, around the complete section. Both results 
are in very close agreement with the datum solutions with a 
tendency for slightly lower pressure inside the duct compared 
with Hess's solution. This could be caused by the slightly 
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Figure 23. Comparison of Calculated y Component of Velocity at 
Two Stations on a Wing with Strake. NACA 0002 Sec- 
tion; a = 5°. 
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REF. 

18 


Figure 24. Comparison of Internal Pressure Calculated using 
a Range of Methods Showing Leakage of Low-Order 
First Generation Model. 


I 


66 



Figure 25. Comparison of Calculated Pressure Distributions on 
a Nacelle with Exit Diameter/Chord .3. NACA 0005 
Section; a = 0°. 
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smaller duct cross-section area generated by the present 10 flat 
panel representation. Again, there is no leakage problem with 
this formulation. 


4.4 W in q~ Fla P Cases 

Cases with slotted flaps have proven troublesome for some 
earlier, low-order panel methods, especially those using a speci 
fied doublet or vorticity singularity distribution and employing 
"Kutta points" where tangential flow is specified just downstream 
of the trailing edges (e.g., (23)). In closely interacting 

multi-element configurations the interference of the downstream 
element acting on the "Kutta point" of the upstream element can 
result in a circulation error; in addition, the solution can be 
sensitive to the shape of the prescribed loading distribution on 
the upstream element — this is quite different from that on the 
regular single-element wing section. 

Two wing-flap cases are considered below for the VSAERO 
program. The first case is one of the exact two-dimensional, 
two-element cases provided by Williams (26). Figure 26 shows the 
close agreement between the VSAERO calculations and the exact 
solution for configuration B on both the main element. Figure 
26(a) and the flap. Figure 26(b). Part of the high aspect ratio 
paneling in the VSAERO calculation is shown in Figure 26(c), and 
has 110 panels chordwise on the main and 60 on the flap. 

The second wing-flap case involves a part-span flap. The 
inset in Figure 27 shows a general view of a wing-flap configura- 
tion at 16 ° incidence. The wake shape is shown after one itera- 
tion cycle. The wing has an aspect ratio of 6 and the basic 
section is a NACA 0012. The 30% chord Fowler flap is deflected 
30° and the flap cutout extends to 0 = .59 on the semispan. The 
surfaces of the wing tip edge, the flap edge and the end of the 
cove are covered with panels. 


Vortex separations were prescribed along the wing tip edge 
and the flap edge. The present calculations of surface pressures 
are compared with experimental data (27) at two stations in 
Figure 27. A section through the flap at n = .33 is shown in 
Figure 27(a) and one just outboard of the flap edge at n - .6 in 
Figure 27 (b) . 


Boundary layer calculations have not been included in this 
calculation yet, and so the agreement between the results at this 
incidence (16°) is probably too close. The lower peak suction in 
the flap data. Figure 27(a), may be caused by interference from a 
nearby flap support wake. 




FLAP: 60 PANELS 


o 



Fgiure 26. Continued. 
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Figure 26. Concluded. 
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Figure 27. Comparison of Pressure Distributions on an Aspect Ratio 6 Wing with 
Part-Span Fowler Flap; a - 16°. 


PRESENT CALCULATION 
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(b) Section at 


4.5 Body Alone 

The case of 4:1 spheroid at 5° angle of attack was run to 
test the prediction accuracy of body cases. Figure 28 shows the 
very good agreement between the VSAERO prediction and the exact 
(non-lifting) solution. The calculation has 10 panels around the 
half-body and 20 longitudinally/ with the case being run with the 
vertical plane of symmetry. Small discrepancies in the C p ~ Z 
plot are probably partly due to the control points being on the 
inscribed polygon (i.e., a plot of Cp 6 would look better, but 
the C p Z plot is more practical for the general case. 


4.6 Wina-Body Cases 

Concave corners have posed a problem to some low-order 
methods in the past, making it necessary to use high panel densi- 
ties. Also, such a region might cause a numerical differentia- 
tion scheme to misbehave. Two wing-body cases are considered 
here. Figure 29(a) shows a longitudinal cut through a pressure 
distribution on a wing-body case from Reference 28 for a = 4° . 
The cut passes just above the wing-body junction. The present 
calculations compare well with other solutions and are in good 
agreement with experimental results. 

A spanwise cut through the present calculations at X/L = .48 
is shown in Figure 29(b) together with the surface geometry 
shape. The pressure distribution passes smoothly from the body 
surface onto the wing surface even though relatively large panels 
are present on the fuselage side (only 6 around the half sec- 
tion) . 

In this calculation the wing and body form one complete 

volume in which $ . = <J> . 

1 00 

The second wing-body case is the Swearingen Metroliner wind 
tunnel model. Figure 30 shows the VSAERO paneling and it will be 
noticed that there has been no treatment of the intersection 
region; i.e., the wing and fuselage panel models are from 
separate wing-alone and body-alone runs. This is not a recom- 
mended practice; however, it is interesting to observe that — at 
least in the present configuration — the results are good and 
compare very well with experimental data right into the junction 
region. (Note that the doublet gradient procedure has not used 
panel information that falls inside the fuselage or wing.) 

Comparisons between experimental (29) and calculated pres- 
sure distributions are made in Figure 31, 32 and 33 at two span- 
wise stations, y = 9.0 (y/b/2 = .22) and y = 18.0 (y/b/2 = .43). 
The location of the body side is at y = 5.0 (y/b/2 = .12). The 
calculations include one wake relaxation iteration. Superimposed 
on each plot is the geometry of the local cross section. 


74 





o 



0) 

a> 

p 

g 

c/) 

05 

a) 

p 

a, 

c 

o 

-P 

V 

G 

G 

>1 

0 
PQ 

1 

tn 

G 

•rl 

s 

a) 

> 

o 

A 

< 


rd 


II 

a 

-M 

rd 

c 

o 

■H 
-P 
rd 
p 
G 
C n 

-H 

LW 

G 

O 

o 

>i 

0 
CQ 

1 

tn 

G 

-rH 

£ 

rd 

G 

O 

05 

G 

0 

•H 

-P 

G 

■H 

P 

■P 

05 

*H 

Q 

CD • 

P 

G • 

05 

05 II 
0) 

P 8 
P-i s 


CM 

<D 

P 

G 

tn 

'H 

(P 


76 


I 


CROSS SECTION SHAPE 


C p VALUES ( 

(PRESENT l 
CALCULA- 
TIONS) 




(b) Spanwise Cut through Present Calculations at 
X/L = .48 

Figure 29. Concluded. 
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Figure 30. Swearingen Metro W.T. Model. 



x/c 

(a) Buttline Cut Y = 9.0. 



«• «8 « 46 48 tl 88 84 

X 



X/C 


(b) Buttline Cut Y = 18.0. 

Figure 31. Comparison of Calculated and Experimental Pressure 

Distributions on the Swearingen Metroliner Wind Tunnel 
Model at a = 0.0, = 2.0. 
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x/c 


(a) Buttline Cut, Y = 9.0. 



(b) Buttline Cut, Y = 18.0. 

Figure 32. Comparison of Calculated and Experimental Pressure 

Distributions on the Swearingen Metroliner Wind Tunnel 
Model at a = 12°, M_ = 2.0. 



(a) Buttline Cut, Y = 9.0. 




Figure 31 shows the comparisons at a = 0°. The experimental 
results at the inboard station indicate a small separation zone 
near the trailing edge and show slightly lower Cp values near the 
leading edge compared with the calculated data. The latter does 
not include a viscous/potential iteration for this case. The 
outboard station shows a much closer agreement between the re- 
sults (Figure 31(b)). 

Figure 32 shows the comparisons at ot = 12 . The results are 
in remarkably good agreement except for the separated zone r 
the trailing edge, which has grown to approximately 20% ot tne 
chord. 

Figure 33 shows the comparison at a = 16°. The calculations 
include one viscous potential iteration and the figure indicates 
the calculated location of separation which is in reasonable 
agreement with experimental data. Figure 34 shows the calculated 
surface streamlines and separation zone on the wing for this 
case. 
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UPPER SURFACE STREAMLINES 
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5.0 CONCLUSIONS 


The formulation of the VSAERO program is described and a 
number of basic test cases examined. The panel method module of 
the program, based on simple quadrilateral panels with constant 
doublet and source, has produced results of comparable accuracy 
to those from higher-order methods given the same density of 
control points. The calculated results compare very well with 
exact, nominally exact and also experimental data cases. Prob 
lems associated with earlier low-order panel methods do not 
appear with the present formulation based on internal Dinchlet 
boundary conditions. Moreover, the low computing cost of the 
simple panel method makes it practical for applications to non^ 
linear problems requiring iterative solutions, e.g., for wake 
shape and surface boundary layer effects. 
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APPENDIX A 


SUBROUTINE FLOW CHART 
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Figure A 1. Flow Chart of Main Program (Deck MVP); 

The Vi scous /Potent ial Iteration Loop. 
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Figure A- 2 . Flow Chart of Subroutine VSAERO in Potential Flow, 
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Figure A-2. Continued. 
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Figure A-2. Continued. 
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Figure A-2. Continued. 
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Figure A-2. Concluded. 




Figure A-3. 


Flow Chart for Major Subroutines in the 
Boundary Layer Analysis (Subroutine BLCTRL) . 
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APPENDIX B 


BIQUADRATIC INTERPOLATION 
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The biquadratic interpolation scheme described below is 
applied in a number of routines in the VSAERO code. Its simple 
multiplier form is very convenient to use and yet it is a "con- 
strained" cubic; i.e., it cannot oscillate wildly. Experience 
with the routine over a number of years in the code has shown it 
to be a reliable method. 

Given a set of position vectors, Pn, n = l,2 ..., N, defining 
a smooth space curve, additional values are interpolated in, say, 
the interval between P n and P n+ i , Figure B-l. The integrated 

on !?„2 ur • en 9 th t s nr to each point from the beginning of the 
curve, i.e., from P 1 , is generated first. For convenience the 
straight segment length is used across each interval; i.e., s = 
P n+l“ p n 1 ' but arc lengths could easily be substituted. n 

Next, two quadratic curves, q^{a,a^ ) passing through points 

n 

P n-1' P n ' P n+1 ' and q 2 ( a f a 2 } passing through points P , P ,, 

n n n+1 

and p n +2' are formed. Figure B-l. 


The normalized interpolation parameter, a, ranges from 0 to 
1 in the n interval, and has value 


where 


a l at P n-1 and a 2 at P n+2 ' 



(3 


n-1 


' s n )/(s n + l 


s ) ; and 
n 


a 2 “ (s n+2 " S n )/ ^ (S n+l ~ S n ) * 


th . 


In the n interval, a linear combination of q. and q 
defines the biquadratic interpolation curve. 1 2 


P(a) - 0^2(0,012 ) + (1 - a)q^ (a, a. ) 
n n 


The biquadratic is, therefore, a cubic, but it is constrained to 
lie between the two quadratic curves. Clearly, the value of a 

for 3 point distance, s, from the start of the curve (but located 
m the nth interval) is 


a 


s/ (s 


n+1 


s n } 
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The form of the interpolation curve can be expressed in 
terms of biquadratic multipliers, Gl, G2, applied to the four 
local position vectors: 

P(a) = P . Gl(a,a, ) + P G2(a,a, , a_ ) 

n-1 1 n 12 

n n n 

+ P n+1 G2(l - a, 1 - a 2 , 1 - a, ) 

n n 

+ * n+2 G1(1 “ >• 

n 


The forms of Gl, G2 are: 


Gl (a,b) = a (1 - a) 2 /{b(l - b) } 

G2(a,b,c) = (1 - a) {1 - a(l - a)/b - a 2 /c) 


These multipliers, based on the linear combination of two 
quadratics, give continuous slope and a piecewise linear — but not 
necessarily continuous—var iation of second derivative across 
each interval. 


The G multipliers can be differentiated to give the tangent 
vector : 

t(a) * P n-1 Hi ( a, a.) + P n H2(a,a 1 ^, a^) 

- P . . H2 ( 1 - a, 1 - a, , 1 - a x ) 

n n 


- p Hl(l - a, 1 - a 2 1 

n+2 n 


where the tangent multipliers are: 

HI (a , b) = (1 - 4a + 3a 2 )/(b(l - b)} 

H2 (a, b ,c) = (4a - 3a 2 - l)/b (3a 2 - 2a) /c - 1 


Note: this is not a unit vector: lt(a)l = 


ds 

da* 



The G multipliers can also be integrated to give the area 
under the curve between the nth and the (n+l)th points: 

*n = Vl +p n F 2 ( “V l2 n ' 

ds 

+ P F2 ( 1 - a , 1 - a- > + P 2 Fid " « 2 } 

n+1 2 n n+z n 

where the integral multipliers are: 

FI (b) = l/{12b (1 - b)} 


F2 (b,c) = (1 - (1/b + l/c)/6}/2 


This assumes that the value of is constant over the interval. 

da 

In this case, therefore, it would be an advantage to use the arc 
length intervals rather than straight line intervals when calcu- 
lating surface distances. 

The F multipliers are for possible future use in the inte- 
gration of surface pressures to obtain forces and moments. 
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BIQUADRATIC CURVE 
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Figure B-2. Biquadratic Interpolation. 




IWNSA 

Nntorvil /V^onaulCS And 
r<\irr ArimnstMton 

1. Report No. 

NASA CR-4023 


Report Documentation Page 


2. Government Accession No. 


4. Title and Subtitle 

PROGRAM VSAERO THEORY DOCUMENT 


A Computer Program for Calculating Nonlinear 
Aerodynamic Characteristics of Arbitrary 
Configurations 


7. Author(s) 


Brian Maskew 


9. Performing Organization Name and Address 

ANALYTICAL METHODS, INC. 

2133 - 152nd Avenue N.E. 
Redmond, Washington 98052 

12. Sponsoring Agency Name and Address 

NASA Ames Research Center 
Moffett Field, California 94035 


3. Recipient's Catalog No. 

5. Report Date 
September 1987 

6. Performing Organization Code 


8. Performing Organization Report No. 

AMI Report 8416 

10. Work Unit No. 


11. Contract or Grant No. 

NAS2- 1 1945 

13. Type of Report and Period Covered 
Final Report 

July 1984 - May 1987 

14. Sponsoring Agency Code 


— — 1 

15. Supplementary Notes 

Point of Contact: Brian Smith (415) 694-5039 

M.S. 247-1 FTS 464-5039 

NASA Ames Research Center 
Moffett Field, CA 94035 

16. Abstract 

The VSAERO low-order panel method formulation is described for the calculation of 
subsonic aerodynamic characteristics of general configurations. The method is based 
on piecewise constant doublet and source singularities. Two forms of the internal 
Dirichlet boundary condition are discussed and the source distribution is deter- 
mined by the external Neumann boundary condition. A number of basic test cases are 
examined. Calculations are compared with higher-order solutions for a number of 
cases. It is demonstrated that for comparable density of control points where the 
boundary conditions are satisfied, the low-order method gives comparable accuracy 
to the higher-order solutions. It is also shown that problems associated with some 
earlier low-order panel methods, e.g., leakage in internal flows and junctions and 
also poor t railing-edge solutions, do not appear for the present method. Further, 
the application of the Kutta conditions is extremely simple; no extra equation or 
trailing-edge velocity point is required. The method has very low computing costs 
and this has made it practical for application to nonlinear problems requiring 
iterative solutions for wake shape and surface boundary layer effects. 


17. Key Words (Suggested by Author(s)) 

Subsonic Panel Method; Arbitrary 
Configuration; Viscous Flow/Potential 
Flow Coupling; Relaxed Wake Iteration 


18. Distribution Statement 


Subject Category 02 


19. Security Classif. (of this report) 

20. Security Classif. (of this page) 

21. No. of pages 

22 . Price 

Unclassified 

Unclassified 

100 



NASA FORM 1626 OCT 86 


NASA-Langley, 1987 



